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Diffusing-Wave Spectroscopy (DWS) extends dynamic light scattering measurements to samples 
with strong multiple scattering. DWS treats the transport of photons through turbid samples as a 
diffusion process, thereby making it possible to extract the dynamics of scatterers from measured 
correlation functions. The analysis of DWS data requires knowledge of the path length distribution 
of photons traveling through the sample. While for flat sample cells this path length distribution can 
be readily calculated and expressed in analytical form, no such expression is available for cylindrical 
sample cells. DWS measurements have therefore typically relied on dedicated setups that use flat 
sample cells. Here we show how DWS measurements, in particular DWS-based microrheology 
measurements, can be performed in standard dynamic light scattering setups that use cylindrical 
sample cells. To do so we perform simple random walk simulations which yield numerical predictions 
of the path length distribution as a function of both the transport mean free path and the detection 
angle. This information is used in experiments to extract the mean-square displacement of tracer 
particles in the material, as well as the corresponding frequency-dependent viscoelastic response. 

An important advantage of our approach is that by performing measurements at different detection 
angles, the average path length through the sample can be varied. Using measurements on a 
single sample cell, this gives access to a wider range of length and time scales than obtained in a 
conventional DWS setup. Such angle-dependent measurements also offer an important consistency 
check, as for all detection angles the DWS analysis should yield the same tracer dynamics, even 
though the respective path length distributions are very different. We validate our approach by 
performing measurements both on aqueous suspensions of tracer particles and on solid-like gelatin 
samples, for which we find our DWS-based microrheology data to be in very good agreement with 
rheological measurements performed on the same samples. 

PACS numbers: 47.57.-s, 46.35.-l-z, 78.35.-|-c, 07.60.-j 


I. INTRODUCTION 


Since its development in the 1980’s, Diffusing-Wave 
Spectroscopy (DWS) [IK3] has proven to be an important 
and versatile tool for studying the dynamics, mechanics 
and structure of a wide range of soft materials. HHII] 
By taking advantage of the fact that the transport of 
photons through an optically turbid sample can be de¬ 
scribed as a diffusion process, DWS extends Dynamic 
Light Scattering (DLS) measurements to the highly mul¬ 
tiple scattering regime. It thus enables access to the 
dynamics of a material at very short time and length 
scales. The method is particularly useful when combined 
with the concept of microrheology, where information on 
the dynamics of tracer particles added to a material are 
used to extract information on the material’s viscoelas¬ 
tic properties. HMH However, the proper analysis of 
any DWS measurement requires detailed knowledge of 
the path length distribution P{s) for photons traveling 
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through the sample to the detector. For a number of sam¬ 
ple geometries and experimental situations, the calcula¬ 
tion or estimation of P{s) has been described in previous 
studies, including for the situation of backscattering from 
a flat sample cell of infinite thickness, or for transmission 
through cone-plate cells or flat circular cells of finite di¬ 
ameter and thickness. [giniiii! Importantly, for sample 
cells in the shape of a flat slab of thickness L, infinitely 
extended in height and width, P{s) can be expressed in 
analytical form, and the analysis of DWS data is there¬ 
fore straightforward. [iin] 

For the cylindrical sample cells used in conventional dy¬ 
namic light scattering setups, however, an analytical ex¬ 
pression for P{s) is not available. DWS measurements 
are therefore usually performed in dedicated instruments 
that use flat sample cells. 

In this paper, we show how DWS-measurements can be 
performed in a standard dynamic light scattering setup, 
using cylindrical sample cells. We perform simple numer¬ 
ical random walk simulations to account for the propa¬ 
gation of photons through a cylindrical cell and describe 
how this information is used to obtain the mean-square 
displacement of tracer particles from the temporal au¬ 
tocorrelation functions determined in experiments. We 
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further show that by performing measurements at dif¬ 
ferent detection angles, the range of accessible time and 
length scales can be extended; this is in analogy to stan¬ 
dard DWS measurements employing a range of different 
cell thicknesses. Importantly, our approach also provides 
a valuable consistency check, especially in the context 
of microrheology, since measurements taken at different 
detection angles should yield the same viscoelastic re¬ 
sponse, even though the corresponding correlation func¬ 
tions must be very different due to the variation in ge¬ 
ometry and average path length. 

In analogy to a conventional DWS measurement, the 
transport mean free path I* is determined by a calibration 
measurement on tracer particles of well known uniform 
size, suspended in a Newtonian liquid of known viscos¬ 
ity; the expected single particle dynamics is thus known 
a priori. Alternatively, the transport mean free path can 
also be directly determined from the measured scatter¬ 
ing intensity as a function of angle, I{9) if an initial 
experimental calibration is combined with results from 
our simple numerical calculations. In the highly multiple 
scattering limit, and in the absence of absorption in the 
sample, I{9) is well approximated by a function that de¬ 
pends only on I*, on the corresponding calculated path 
length distribution P{s), and on a constant /3exp that is 
determined by the experimental setup. 

We thus demonstrate that standard goniometer light 
scattering setups can be successfully used for DWS mea¬ 
surements. Compared to dedicated DWS setups, our 
method has the advantage of being able to reliably deter¬ 
mine the transport mean free path I* as well as to extend 
the range of accessible length and time scales, using only 
a single cylindrical sample cell. 

We illustrate and test the use of our approach by per¬ 
forming DWS-based microrheology measurements on a 
typical solid-like soft material, gelatin, and find the re¬ 
sulting frequency-dependent viscoelastic moduli to be in 
very good agreement with separately performed rheolog¬ 
ical measurements. 


II. EXPERIMENTAL: MATERIALS AND 
METHODS 

A. Background on DLS, DWS, and microrheology 

Standard static and dynamic light scattering experi¬ 
ments are limited to samples that exhibit very little mul¬ 
tiple scattering, with the overwhelming majority of de¬ 
tected photons having been scattered only a single time 
within the sample. A typical setup for such single scat¬ 
tering experiments uses a cylindrical sample cell that is il¬ 
luminated by a laser, as shown schematically in FigQA). 
The detector, typically comprising an optical fiber that is 
coupled to a photomultiplier tube, can be positioned at a 
range of detection angles 0, corresponding to scattering 
wave vectors q{9) = sin(0/2), where n is the refrac¬ 

tive index of the sample and A is the wavelength of the 



FIG. 1: (Color online) Schematic experimental setup for typ¬ 
ical dynamic light scattering measurements. (A) Standard 
dynamic light scattering (DLS) setup employing a cylindrical 
sample cell and a goniometer, which enables accessing dif¬ 
ferent scattering angles 0. (B) Diffusing-Wave Spectroscopy 
(DWS) setup in transmission geometry. The pathways of pho¬ 
tons are well-described by a random walk of step size I*. (C) 
Schematic of random walk simulations in a cylindrical geome¬ 
try. The detection angle 9 is defined as for conventional DLS; 
here it determines the distance L{9) between the points of en¬ 
try and exit of detected photons. (D) Expected distribution 
p{z) of the a coordinate where a photon exits the cylinder; a 
fraction cxz (marked area) reaches the detector. 


laser in vacuum. For single scattering, the fluctuations in 
the detected intensity, which reflect the dynamics of the 
scatterers, are then quantified by the temporal intensity 
autocorrelation function 
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where t is the lag time and the brackets < ..>i indicates 
a time-average over all times t. The field autocorrelation 
function gi(t), measured at a wave vector g, reflects the 
temporal fluctuations of the electric field. It can be re¬ 
lated, via the so-called Siegert r elation to the intensity 
correlation function, as g\(t) « a/( 52 (t) — 1) //?, where (3 
is the coherence factor, m a constant that depends on 
the experimental setup. For a Gaussian distribution of 
displacements Ar, the field correlation function gi{t) is 
directly linked to the dynamics of scatterers in the sam¬ 
ple, as 


gi(t)=e-^(^'-'(‘)), (2) 

where (Ar^(t)) is the time-dependent mean-square 
displacement of scatterers in the material. For the sim¬ 
plest example, where the scatterers are uniformly sized 
particles suspended in a Newtonian liquid, the particles 
undergo ideal Brownian motion, and thus (Ar^(t)) = 
6Dt, where D is the particle diffusion coefficient in 3 di¬ 
mensions. For this case, the field correlation function 
has a single exponential form, gi{t) = e”'"*, where the 
q-dependent decay rate F = Dq^ is set by D. 

Dijfusing-Wave Spectroscopy (DWS) is an extension 
of dynamic light scattering measurements to the highly 
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multiple scattering regime. A typical experimental setup 
for DWS is shown in Figj^B). In contrast to conven¬ 
tional DLS measurements, this technique requires that 
photons are scattered many times, before they reach the 
detector. For this highly multiple scattering regime, the 
propagation of photons through the sample can be ade¬ 
quately described as a simple diffusion process, where the 
details of each single scattering event are no longer rele¬ 
vant. This photon diffusion processes can be accounted 
for by a single parameter, the so-called transport mean 
free path I*. This characteristic length scale is defined as 
the average distance a photon travels in the sample be¬ 
fore its direction of propagation is randomized. The path 
of photons through the sample can thus be approximated 
as an ideal random walk with step size I*. For such a ran¬ 
dom walk the path length of photons, and the number 
of randomizing scattering events, is no longer uniform, 
as is the case for single scattering. Instead, the correla¬ 
tion function measured in an experiment is determined 
by contributions from all path lengths s weighted by the 
path length distribution P{s), as 

Z*®® u 2 

g,{t) = / + (3) 
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where fcg = 2TTnjX is the magnitude of the photon wave 
vector in the sample and s/l* reflects the number of ran¬ 
domizing scattering events experienced by a photon with 
path length s. [3] The basis for this simple form of Eq|^is 
that each of the approximately s/l* randomizing scatter¬ 
ing events contributes to a change of this particular pho¬ 
ton path by a squared distance of (Ar^(t)), leading to 
a partial decorrelation of gi{t). The cumulative decorre¬ 
lations from all these randomizing scattering events thus 
lead to the functional form in Eq|^ Knowledge on the 
path length distribution P{s) is therefore essential in the 
analysis of DWS measurements; without such knowledge 
the measured correlation functions cannot be related to 
the dynamics of the scatterers. The path length distri¬ 
bution depends sensitively on the geometry of the sam¬ 
ple cell used in the experiment. For sample cells in the 
shape of a flat slab, infinitely extended in both height and 
width, P(s) can be expressed in analytical form as a func¬ 
tion of I* and the thickness L of the sample cell. lain] 
This is one of the main reasons why DWS measurements 
have typically relied on measurements performed in ded¬ 
icated instruments, employing flat sample cells. 

Such dedicated DWS instruments can also offer other 
important advantages, in particular for measurements on 
solid-like, nonergodic samples, where the measured, time- 
averaged correlation functions are not representative of 
the ensemble-averaged dynamics of the sample. [20] 
Methods for acquiring ensemble-averaged correlation 
functions in DWS measurements include the use of 
double-cell techniques, where either an ergodic sample 
with slow dynamics [21] or a slowly rotating opaque 
disc [12] is placed in front of the sample cell. Both these 
techniques create a slow randomization of the incoming 


photon paths, resulting in an ensemble-averaging of the 
collected temporal correlation functions. Either transla¬ 
tions of the sample cell or rotations of an opaque disc 
can also be employed for ensemble-averaging using echo 
techniques, jUj yielding ensemble-averaged correlation 
functions at long time scales and with excellent statis¬ 
tics. HU HU While these ensemble-averaging techniques 
could in principle also be incorporated into a standard 
goniometer setup, we choose an alternative method, so- 
called Pusey-averaging. This method uses the measured 
time-averaged correlation function and the measured 
ensemble-averaged scattered intensity together with a 
simple theoretical treatment to provide the ensemble- 
averaged correlation function. |25l |26] 

The ensemble averaged scattering intensity (J)^ can be 
readily acquired in separate intensity measurements dur¬ 
ing which the sample is rotated; the measured dynamics 
is perturbed by the motion of the sample, but the average 
scattering intensity is still properly ensemble-averaged. 
Using the ratio of the ensemble-averaged to the time- 
averaged scattering intensities Y = the ensemble- 
averaged field autocorrelation function gi (t) can then be 
estimated as a function of the time-averaged correlation 
function as 

= + (4) 

where 52 (t) = 1 + is the time-averaged inten¬ 

sity autocorrelation function normalized by the coherence 
factor /I, which is obtained from the separate ensemble- 
averaged measurements, and cr^ = ^ 2 ( 0 ) — 1 characterizes 
the short-time intercept of 52 (t) • 

We can now use the resulting ensemble-averaged field 
autocorrelation function gi{t) to extract viscoelastic 
properties of the sample, using the microrheology con¬ 
cept mm- To do so, we employ the local power-law 
approximation developed by Mason et al jUKle]. In 
brief, the method is based on the assumption that the 
Stokes-Einstein relation, which links the thermal motion 
of particles in a Newtonian liquid to the viscosity of the 
surrounding liquid, can be generalized to viscoelastic ma¬ 
terials with frequency-dependent linear viscoelastic mod¬ 
uli. The approximation also neglects inertial effects on 
the motion of the probe particles, which is justified for 
most soft materials at frequencies below « I MHz. 

By describing the time-dependent mean square dis¬ 
placement as a local power-law around each data point, 
the magnitude of the frequency-dependent complex mod¬ 
ulus can be expressed in analytical form as 

knT 

7 ra(Ar2(l/w))r(l + a(I/w))’ 

where a is the particle radius, k^T the thermal en¬ 
ergy, a{t) = is the logarithmic slope of the 

mean-square displacement as a function of lag time, and 
F denotes the gamma function. 






4 


B. Sample preparation 

Polystyrene particles (micromod Partikeltechnologie 
GmbH, Germany) coated with a grafted layer {M^ = 
300 g/mol) of poly (ethylene glycol) were used as tracer 
particles in the DWS measurements. The diameter of 
the particles is 1 /rm and they are provided suspended 
in water at a concentration of 5 wt%. The test samples 
with tracers in water are prepared by mixing the stock 
particle suspension with deionized water (Milli-Q water, 
a > ISMfl • cm at 25 G), to obtain the desired tracer 
particle concentrations Ctracer- To study the effect of the 
transport mean free path G, which is expected to scale 
as I* oc 1/Ctracer, we prepare a series of samples with 
particle concentrations ranging from Ctracer ~ 0.3 wt% 
to 5 wt%. 

The aqueous gelatin gel is prepared by mixing water 
with 5 wt% gelatin powder (type A, from porcine skin, 
Sigma, USA) and 1.25 wt% of tracer particles at elevated 
temperatures of « 60° G. The mixture is homogenized 
for 30 minutes using a magnetic stirrer, transfered to 
the cylindrical sample cell used in the experiment, and 
subsequently allowed to cool down to room temperature. 

C. Light scattering experiments 

All dynamic light scattering experiments are per¬ 
formed in a static and dynamic light scattering setup 
(ALV GGS-3, ALV GmbH, Germany), equipped with a 
50 mW solid state laser (A = 532 nm) and a goniome¬ 
ter that allows for variation of the detection angle from 
0 « 20 deg to 0 « 160 deg. Measurements are performed 
in cylindrical cells with outer diameter 10 mm and inner 
diameter 8.65 mm; the cell radius relevant to the propa¬ 
gation of photons in the sample cell (see FigQC)) is thus 
R « 4.33 mm. Measurements of 30 s duration are per¬ 
formed at detection angles between 30 deg and 150 deg 
in steps of 10 deg. To minimize the detection of stray 
light, reflected from surfaces in the setup, our measure¬ 
ments are performed in vertical-horizontal mode, with 
the incoming light vertically polarized and a horizontal 
polarizing filter placed in front of the detector. 

For the gelatin samples we use separate experiments 
on the same sample to determine the ensemble-averaged 
scattering intensities (1)^ as well as the coherence fac¬ 
tor fj needed for the Pusey averaging method. In these 
separate experiments the sample cells are slowly rotated 
during data acquisition; we perform three such measure¬ 
ments at each scattering angle, each lasting 10 seconds. 

In principle, the light scattering measurements we de¬ 
scribe in this paper could be performed with any stan¬ 
dard goniometer setup; however, not all laser sources 
that are included in standard goniometer setups may be 
suitable for performing DWS measurements. In partic¬ 
ular, as a result of the long path lengths of the pho¬ 
tons through the sample, DWS requires a laser source 
with a sufficiently long coherence length. Information 


on the coherence length of a laser is difficult to obtain 
from the standard information provided by manufactur¬ 
ers, and its measurement requires complex setups only 
available in specialized optics laboratories. Within the 
range of concentrations studied here, the longest relevant 
path lengths of photons through the samples were limited 
to around 3 meters. This would be a problem for instance 
if a Helium-Neon laser were used, which has a typical co¬ 
herence length of only around 20 cm. Semiconductor 
lasers, and lasers coupled into single mode fibers, how¬ 
ever, typically have much longer coherence lengths that 
can reach hundreds of meters. While we have not directly 
measured the coherence length of our laser source, the 
good agreement of our path length simulations with ex¬ 
periments (see Results section) makes us confident that in 
our setup the coherence length of the laser is larger than 
the relevant path lengths of photons traveling through 
the sample. 


III. RESULTS AND DISCUSSION 
A. Simulation of photon paths through the sample 

To properly interpret experimental data in a setup 
with a cylindrical cell, the path length distribution P{s) 
of photons traveling through the sample is required both 
as a function of the detection angle 6 and the transport 
mean free path I*. 

To achieve this, we perform numerical simulations of 
photons traveling through a cylindrical cell, assuming 
that they undergo an ideal random walk with step size I*. 
In the 2-dimensional coordinate system given in FigQC), 
photons are released at point {x/R = —1 -I- I* /R,'^R = 
0), where I* is the transport mean free path and R is the 
radius of the cylindrical cell. Subsequently, each photon 
is propagated in steps of U/i?, where each step proceeds 
in a random (3D)-direction. At the point where the pho¬ 
ton exits the cell ( ), we evaluate the number 

M of scattering events, and record the detection event 
with respect to the observed detection angle 9. 

We do this by dividing the surface of the cylindrical 
cell into Ubins angular bins, spanning from 0 to 180 deg 
(taking into account the symmetry around the x-axis). 
In addition, to take into account the 3-dimensional na¬ 
ture of photon transport in the real geometry, we consider 
the fact that each realization of a 2-dimensional photon 
path represents a whole range of possible 3-dimensional 
paths with an identical number of scattering events M 
and identical (x, y)-paths. Since in the ^-direction the 
photon also performs a (1-dimensional) random walk, we 
can readily express the probability distribution p{z) for 
the photon to end up at a position z after propagating 
M random steps. What is relevant here is the fraction 
of those paths that will reach the detector, as illustrated 
in FigQD). We assume that all photons with \z\ < Az 
are detected; in accord with the resolution of the angu¬ 
lar bins, we set Az = —. Then, the fraction of 
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FIG. 2: (Color online) Simulation results. (A) Path length distribution P{s) for different values of I* /R, calculated at a fixed 
detection angle 0 = 90 deg. (B) Master curve of scaled path length distributions, showing that the shapes of P(s) calculated 
at different I*/R are similar. (C) Corresponding average path lengths s as a function of R/R. As expected, we find a scaling 
s oc 1/Z*, as indicated by a power-law fit to the data (dashed line), yielding an exponent m = —0.975 ± 0.05. (D) Average 
number of scattering events s/l* as a function of the distance L{6) between the entry and exit points of the detected photons. 
The dashed lines serve as a visual reference, indicating the scaling s/R oc L(0)^ that would be expected for an unrestricted 
random walk. 


contributing 3-dimensional paths is given as 

u ^ fa\ 

where erf is the error function, and is the effec- 

tive 1-dimensional step size in the z-direction. We thus 
account for diffusion in the z-direction in our statistics of 
path length distributions by, instead of adding 1, adding 
a contribution to the angular bin corresponding to 
each simulated photon path: /(ribin) —t finun) + cTz- 
Each bin thus represents a detection area of surface 

area Abin ~ cumulative value of each angu¬ 


lar bin, after propagating N photons and normalizing by 
fV, thus defines a (dimensionless) scattering intensity as 
dsim := , representing the probability for a photon 

to reach the detection area corresponding to bin number 

^bin- 

In addition to recording the angle where the photons 
end up, we also record, for each angle 0(nbin), a distri¬ 
bution of the number of scattering events, by adding a 
contribution to a bin accounting for the number of 
scattering events at each angle 6*(ribin)- The bins are lin¬ 
early spaced, with bin number 100 representing a num¬ 
ber of approximately {L{9)/R)'^ scattering events, which 
corresponds to an expected average number of scatter¬ 
ing events for a distance L{9) between the entry point 
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and the detection point of the photons.^ We use 300 
bins per angle d(nbin), thus accounting for up to 3 times 
the expected typical number of scattering events; higher 
numbers, while not counted, in practice are extremely 
rare in our simulations and do not significantly affect the 
resulting path length distributions. 

In order to achieve good statistics in the calculated 
path length distributions, the paths for a large number 
of photons have to be simulated. 

In the actual experiments obtaining good statistics is 
usually not a problem, due to the enormous number of 
photons that are propagated. In our experiments we use 
a laser with 50 mW of power at a wave length of 532 nm; 
this corresponds to « 10^^ photons entering the sample 
cell, every second. Such numbers are beyond the capabil¬ 
ity of computer simulations; in comparison, for our cal¬ 
culations we typically simulate 10® photons propagating 
through the sample, which is enough to yield reasonable 
statistics, and relatively smooth calculated path length 
distributions. 


B. Scaling properties of P{s) 

Typical obtained simulation results for P(s) are shown 
in Figj^A); these curves are calculated at a fixed angle 
0 = 90 deg for different values of l*/R. Interestingly, 
while the average path length decreases with increasing 
I*, the shapes of these path length distributions appear 
surprisingly similar, . 

In fact, we can overlay the curves from Fig. [^A) and 
create a master curve, as shown in Fig. [^B). To obtain 
this master curve, we have rescaled the path length with 
a factor s and multiplied the magnitude with the same 
factor; it turns out that s is the average path length, 
defined below in Eq. 

Any practical use of the calculated path length dis¬ 
tributions requires that P{s) data are available for any 
arbitrary value of I*. To address this problem, we calcu¬ 
late path length distributions for different values of I*/R 
and examine the scaling properties of these path length 
distributions. 

In contrast to an ideal random walk, the path of pho¬ 
tons through our cylindrical sample is constricted by the 
geometry. Nevertheless, the essential scaling properties 
of a random walk still hold approximately for the path 
length distributions simulated here. In particular, for 
an unrestricted random walk of step size /*, we expect 
the mean square displacement (^AR'^') to be given as 
(Ai?®) = Ml*"^, with M the number of steps. Consider¬ 
ing the average path length 


^ In the cylindrical cells studied here, the average pathlengths s 
depend on the detection angle and are typically shorter than 
estimated from s/l* ^ as seen in Fig&B). 


s:= sP{s)ds , (7) 

Jo 

we can estimate the average number of scattering 
events M to travel to a point at distance L(6) from the 
origin to be approximately given as M ~ {L{0)/1*) , as 
would be the case for a completely unrestricted random 
walk. 

As the path length is s = Ml*, the average path length 
should scale with L{9) and Z* as s « L{9)'^/1*. Con¬ 
versely, at fixed detection angle 9 and thus fixed distance 
L{9), we would clearly expect a scaling of s cx l/l*. 

To test this scaling, we examine the P{s) data with 
respect to both I* and the detection angle 9, where a 
variation of the latter corresponds to a variation of the 
distance L{9) between the entry and detection points of 
the photons. In FigMC) we plot the average path length 
s as a function of 1*JR, for simulation data calculated at 
a single detection angle 0 = 90 deg. Indeed, the data 
is in excellent agreement with a scaling of s oc 
the dashed line in Figj^C) shows a power-law fit to the 
data, yielding an exponent of —0.975 ±0.05. This scaling 
is a consequence of the self-similarity of random walks, 
which enables us to approximate each random walk with 
a “coarse grained” version of larger step size; this scaled 
random walk essentially follows the same path, but, as a 
result of the increased step size, exhibits a reduced con¬ 
tour length. 

In contrast to this simple scaling as a function of I*, 
if we examine the average number of scattering events 
as a function of L{9), we find significant deviations from 
the naively expected scaling s/l* oc L(9)^, as shown in 
Figj^D). The symbols in this figure show the simulation 
data for different values of fixed I*/R, and the solid lines 
show the simple prediction discussed above, a power-law 
with exponent 2. In hindsight, it is clear that such de¬ 
viations should be expected, as, in contrast to the l*- 
dependence at fixed detection angle, a variation of L{9) 
implies a significant modification of the effective sample 
geometry. It is thus evident that calculations of P(s) 
at different detection angles are necessary. However, the 
simple scaling properties with respect to I*, highlighted in 
FiggB), can be exploited to obtain accurate path length 
distributions for arbitrary Z*-values, based on simulations 
performed at a single value of I*/R. 

C. Determination of I* for samples with known 
tracer dynamics 

Typically, when DWS experiments are used to perform 
microrheology measurements, uniformly sized tracer par¬ 
ticles are added to the soft material of interest. If the 
transport mean free path I* is known, the dynamics of 
these particles can be extracted, yielding direct infor¬ 
mation on the viscoelastic properties of the surrounding 
material. Vice versa, if the viscoelastic properties of the 
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FIG. 3: (Color online) Simulation results versus experiments for colloidal particles in water {(j) = 0.625%, d — Ipm): Measured 
(blue circles) and calculated (red line) field correlation function gi{t) for a detection angle of 0 = 50 deg (A), 6^ = 90 deg (B), 
and d — 130 deg (C). ln(gi(t)) as a function of t for the same values of 6 = 50 deg (D), 6 — 90 deg (E), and 6 — 130 deg 
(F). (G): I* obtained from fitting simulation data to experiments, as a function of 6 for different (f>. (H): The same I* values 
plotted as a function of (p for different values of the detection angle 9-, we find approximately I* oc !/<)>, as illustrated by the 
dotted line. The data is also in fair agreement with predictions from Mie scattering theory. [271429] shown as a dashed line. 
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FIG. 4: (Color online) Obtaining Z* from intensity measurements: Calculated and measured intensities as a function of 9 and 
I*. (A) Intensity as a function of 9 for simulations employing different values of I*. (B) Concentration dependence of the 
measured scattering intensity (for different angles 9). (C) Comparison between experimental and simulation results; scaled by 
a single constant /3exp good agreement between simulation and theory is obtained at all concentrations studied. Thus, once 
/lexp is determined, the measured transport mean free path Z* can be directly deduced from 1{9). 

























surrounding material are a priori known, then I* can be 
extracted from a DWS measurement. A requirement is 
that the particles scatter much more strongly than the 
material, ensuring that any detected dynamics are re¬ 
lated only to the particle dynamics, and not to fluctu¬ 
ations within the surrounding sample. If this criterion 
is fulfilled, the simplest method for determining I* is to 
measure the scattering of uniformly sized, spherical par¬ 
ticles, suspended in a Newtonian background liquid of 
known viscosity. 

To test our approach, and for calibration of the trans¬ 
port mean free path I*, we here perform a series of ex¬ 
periments using samples with different concentrations of 
uniformly sized polystyrene particles (coated with poly¬ 
ethylene glycol, Mw « 300 g/mol, 1 pro diameter, pur¬ 
chased from micromod GmbH, Germany) suspended in 
water. We measure the field autocorrelation function 
gi{t) of the scattered light for these samples and test 
how the correlation functions predicted from our pho¬ 
ton path simulations compare to these experimental data. 
As shown in Figj^A-G), where we show data on a sus¬ 
pension of particles at a volume fraction </> = 0.625% 
and detection angles 0 = 50 deg, 90 deg, and 130 deg, 
we obtain remarkably good agreement between exper¬ 
iments (shown as blue circles) and simulations (shown 
as red lines), where I* is the only adjustable parame¬ 
ter. While the dynamics is expected to be purely Brow¬ 
nian, with (Ar^(t)) increasing linearly with time t, due 
to the broad path length distribution of photons passing 
through the sample, gi{t) deviates significantly from the 
single exponential decay that would be observed in single 
scattering experiments. This can be more clearly seen 
in Figj^D-F), where ln(( 7 i(t)) is plotted as a function 
of time; in such a plot an exponential decay would ap¬ 
pear as a straight line, as illustrated by the dashed lines, 
which show exponential hts to the short-time regime of 
gi(t). The non-exponential shape of the data thus be¬ 
comes evident and is captured very well by the curves 
predicted from our simulations, shown as red lines. Fits 
performed for the same sample, but at different angles, 
should yield the same G-values. Indeed, we obtain good 
agreement between the G-values extracted from the data 
in Fig[^A-G): we obtain I* = 540 pm, I* = 533 pm, and 
I* = 523 pm at angles of d = 50 deg, 90 deg, and 130 deg, 
respectively. In fact, we obtain good agreement between 
measurements taken at different angles for all the con¬ 
centrations studied, with volume fractions ranging from 
4> = 0.313% to 2.5%. As shown in Figj^G), the fitted 
Z*-values as a function of 0 exhibit only small variations. 
Somewhat larger deviations are observed for the sample 
with the lowest concentration, at the largest detection 
angles. We attribute this to the fact that this sample 
has the longest I*, combined with the shortest distances 
L(0) between entry point and exit point of the photons; 
Z* « 1 mm and L{9) « 2.2 mm and thus L{9)/1* « 2.2. 
In this case the path length of photons is no longer ade¬ 
quately described as an ideal random walk. 

Besides these discrepancies at small values of L{9), the 


fitted Z*-values depend only on the volume fraction, ir¬ 
respective of the detection angle. To examine the (jr- 
dependence of the data, in Figj^H) we plot Z* as a func¬ 
tion of 4>, observing approximately the expected scaling 
I* oc !/(/), 130] as indicated by the dotted line. The data 
is also in fair agreement with Mie scattering calculations 
plotted as a dashed line in Figj^H). [331 HO] The calcu¬ 
lations are performed using the web application available 
on the website of LS Instruments, Switzerland, m using 
as input parameters the particle size, the wavelength of 
the laser A = 532 nm, as well as a refractive index of 
nps = 1.598 for the particles and riHaO = 1-33 for water. 


D. Obtaining I* from intensity measurements 

In the absence of absorption, the intensity detected 
at each angle should be fully determined by the trans¬ 
port mean free path of photons in the sample. While 
absorption is relatively straightforward to include in the 
current data analysis, here we choose to neglect its ef¬ 
fects since absorption is relatively weak in the aqueous 
samples studied; the typical absorption length is much 
longer than the typical path length of photons through 
the samples. As a result, after calibration using a refer¬ 
ence sample of known dynamics, a simple measurement 
of the scattering intensity at different angles on the sam¬ 
ple of interest is sufficient for determining its I*. To vali¬ 
date this, we compare the scattering intensities predicted 
from the simulations with those measured in experiments 
performed on our polystyrene suspensions. 

In Figj^A) we plot the scattering intensity /aim as a 
function of detection angle 9, as predicted from the pho¬ 
ton path simulations. The shape of these curves is very 
different from those typically obtained in single scatter¬ 
ing experiments on dilute suspensions, where generally 
the intensity is highest at small detection angles, corre¬ 
sponding to low q — values. In the highly multiple scat¬ 
tering regime, however, the intensity is generally highest 
for detection points closest to the entry point of photons 
into the sample, which corresponds to large 0-values. 

Importantly, we find that the angular dependence of 
the recorded scattered intensity for the tracer particle 
suspensions agrees remarkably well with the behavior 
predicted from our simulations, as shown in FigQB). For 
both simulations and experiments, we observe a ratio of 
« 40 between the intensities measured at 0 = 30 deg 
and 9 = 150 deg. Moroever, comparing Figj^B) with 
Fig|^A), we observe that the shapes of the simulated in¬ 
tensity curves are very similar to those of the experimen¬ 
tal data. In fact, the two data sets can be superposed 
simply by scaling the simulated curves with one single 
factor /3expj the value of which depends on experimental 
parameters such as the size of the detection area, and 
the distance between the detector and sample. We find 
that good agreement between the measured and simu¬ 
lated curves is obtained for a value of /3exp ~ 8 • 10^° Hz, 
as shown in FigQG). 
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Thus, given /3exp, a good estimate of I* can be de¬ 
termined for a sample of unknown properties solely by 
measuring the scattered intensity at different angles. 


E. Test on a viscoelastic material 

Finally, we test the use of our method in the context 
of microrheology, where the measured correlation func¬ 
tions and the corresponding tracer particle mean-square 
displacements (Ar^(t)) are used for the determining vis¬ 
coelastic properties of a sample. As a test material we 
use a common solid-like soft material, an aqueous gelatin 
gel at a concentration of 5 wt%. The plateau storage 
modulus of this material is on the order of 1 kPa, which 
means that for the case of micron-sized tracer particles 
we need to be able to access particle displacements at 
sub-nanometer length scales. DWS is ideally suited for 
this, since the displacements of all the tracer particles 
encountered by a photon on its path through the sam¬ 
ple cumulatively contribute to changing the total photon 
path length. 

We perform measurements on the gelatin sample for 
detection angles ranging from 0 = 30 deg to 0 = 150 deg. 
The sample is highly non-ergodic, as indicated by the 
fact that the intercept of the measured intensity corre¬ 
lation functions gi{t) varies significantly between mea¬ 
surements. We therefore use the Pusey-averaging proce¬ 
dure to obtain a good estimate of the ensemble averaged 
correlation functions from the measured, time-averaged 
correlation functions, as outlined in the experimental sec¬ 
tion. As expected, and shown in Figj^A), the resulting 
ensemble-averaged field correlation functions gi{t) vary 
with the detection angle 0. These correlation functions 
do not decay significantly; they reach a plateau at values 
of gi{t) > 0.9 at the longest time scales accessed in the 
experiments. This reflects the fact that the gelatin sam¬ 
ple has a relatively high modulus and the thermal motion 
of the tracer particles is therefore limited to short length 
scales. Using the procedure outlined above, we obtain the 
transport mean free path I* of the gelatin sample directly 
from the measured scattered intensities, using the inten¬ 
sity scaling factor /3exp as determined from the measure¬ 
ments on our pure tracer suspensions. The mean-square 
displacements of tracer particles in the gelatin sample 
are obtained by numerically inverting Eqj^ using the 
measured gi{t), Z*, and the calculated 0-dependent path 
length distribution P{s) as input. The corresponding 
mean-square displacements are shown in Figj^B). Given 
the highly nonergodic nature of the sample studied, the 
data taken at different detection angles are in fair agree¬ 
ment; note that the magnitudes of the accessed particle 
displacements are in the sub-nanometer range. We can 
now convert these data to viscoelastic moduli, using the 
microrheology concept lain] and the local power-law 
approximation [HI [IS], developed by Mason et al. The 
magnitude of the resulting complex modulus |G*(a;)| is 
on the order of 1 kPa and depends only weakly on fre¬ 


quency, as shown in Figj^C). The curves obtained for 
different detection angles exhibit significant variations, 
as shown in the inset, where we plot the low frequency 
plateau values G* = \G*{uj = 5 rad/s)| as a function of 
detection angle (the frequency of 5 rad/s is indicated as 
a dotted line in the main plot). Since |G*(a;)| is approx¬ 
imately inversely proportional to the mean square dis¬ 
placement = l/w)), these variations in the mag¬ 

nitude of the complex modulus directly reflect those ob¬ 
served in the tracer mean-square displacements. 

A simple error analysis (see SI) suggests that the main 
sources of errors are on the one hand statistical errors in 
the intensity correlation function as a result of the finite 
measurement duration, and on the other hand errors in¬ 
troduced via the Pusey averaging procedure via the error 
in the intensity ratio Y = between the time-averaged 
and the ensemble-averaged scattering intensities. 

The total corresponding relative error in the modulus, 
AG/G, can be expressed as a function of the relative er¬ 
rors in the intensity ratio, AU/U, and the ratio between 
the probed time scale t and the measurement duration 
T, as [3T| 

AG A(Ar2(f)) AT 3 [T 
G ^ (Ar2(t)) y g,Hnig,)]j T 

We can estimate the error in determining the in¬ 
tensity ratio Y from the standard deviation of the 
3 ensemble-averaged intensity measurements taken at 
each angle as AY « Y ■ A/g/Je, where A/g is taken 
as the standard deviation of the 3 intensity measure¬ 
ments, and Jg is the average ensemble-averaged in¬ 
tensity. The relative error in the mean square dis¬ 
placement cmsd = ^((^^^(^)))/(^^'^(0) that results 
from the error in Y can be estimated as cmsd ~ 
d{Ar^{t))/dgi{t)-dgiit)/dY-AY/Y. As 9gi(t)/5r « 
1/y^ and d{Ar^{t))/dgi{t) « - {Ar^{t))/[ln{gi{t)) ■ 

gi{t)], we find emsd « - [\n{gi{t))gi{t)Y'^] ^ Be¬ 
cause the modulus |G*(a;)| is essentially given as the in¬ 
verse of the mean-square displacement, it has the same 
typical relative error ec* ~ cmsd ■ For our measurements 
we find typical values ec* ~ 0.5, as shown in the inset of 
Figj^C), where the corresponding errors bars are shown 
for each angle. 

While the magnitude of these errors is indeed signif¬ 
icant, the observed angle-dependent variations remain 
within a factor of « 2 and exhibit no clear systematic 
dependence on the detection angle. The latter indicates 
that the observed variations are the result of random, 
rather than systematic measurement errors, and thus 
do not reflect a systematic issue with the data analy¬ 
sis method. This suggests that performing an average of 
the data measured at different detection angles provides 
a more reliable result than data from a single detection 
angle measurement. We thus average the mean-square 
displacements obtained at all accessed angles, using the 
inverse of the estimated error bars as weights in the av¬ 
eraging procedure. Using this averaged data to calculate 
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FIG. 5: Measurements on a solid-like, non-ergodic sample, an aqueous gelatin gel at a concentration of 5 wt% with embedded 
tracer particles (1.25 wt%, a = 1 ). (A): Pusey-averaged field autocorrelation functions gi{t) as a function of lag time t, 

measured for detection angles ranging from 0 = 30 deg to 6 = 150 deg. Curves are offset vertically by increments of -0.05 for 
clarity. (B): Mean-square displacements extracted from the same data, plotted as a function of time, yielding fair agreement 
between measurements taken at different angles. Note that sub-nanometer displacements are accessed. (C) Magnitudes of the 
corresponding complex shear moduli |G*(cu)| as a function of frequency lj. The inset shows plateau values Gp, accessed at a 
frequency of 5 rad/s (indicated as a dotted line in the main plot), as a function of detection angle. Within the estimated error 
bars, no strong trend in the data is observed; averaging over data from different 6 thus appears justified. (D) Storage modulus 
G'{uj) (solid squares) and loss modulus G"(u)) (open squares) averaged over all measurements taken at different angles as a 
function of frequency ui. Comparing these averaged moduli to results from conventional oscillatory rheology, with G'{u}) shown 
as solid black circles and G"{uj) shown as open black circles, we observe very good agreement. 


the viscoelastic response of the sample we obtain an aver¬ 
aged viscoelastic response, shown in Fig. j^D), where we 
plot the storage modulus G' (solid squares) and the loss 
modulus G" (open squares) as a function of frequency w. 
The oscillations observed in the data at a frequency of 
« 50 — 200 Hz are likely the result of a mechanical dis¬ 
turbance or vibration that we could not eliminate in our 
experimental light scattering setup on this highly out-of- 
equilibrium sample. The effect can be directly observed 
in the measured correlation functions at the correspond¬ 
ing time scales, as seen in Fig[^A). For mechanically 
weaker samples we have not observed these types of os¬ 
cillations in our setup. As a result of the importance 


of the time-derivative of in determining the vis¬ 

coelastic moduli, the oscillations in the gi (t)-data are am¬ 
plified in the corresponding viscoelastic moduli. The fre¬ 
quency range where we do not trust the data as a result of 
these oscillations is indicated with a grey background in 
Figj^D). Nevertheless, besides these oscillations, we ob¬ 
tain very good agreement with measurements performed 
using a conventional oscillatory rheometer, shown in the 
same figure as solid black circles for G' and open black 
circles for G". 
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IV. CONCLUSIONS 

We have developed a simple method for properly inter¬ 
preting dynamic light scattering data from highly mul¬ 
tiple scattering samples using a standard dynamic light 
scattering setup with a cylindrical sample geometry. By 
performing ideal random walk simulations within a cylin¬ 
drical geometry, we predict the path length distribution 
P{s) of photons passing through the sample cell. This en¬ 
ables us to extend the use of DWS measurements to stan¬ 
dard dynamic light scattering instruments. The method 
can be applied in the context of microrheology, where 
the dynamics of embedded tracer particles are used to 
access the frequency-dependent viscoelastic response of 
soft materials. 

The main strength of our approach, besides not re¬ 
quiring a dedicated instrument, lies in the fact that by 
varying the detection angle we can access a wide range 
of different effective sample geometries with different av¬ 
erage path lengths, using one single cylindrical sample 
cell. This variation of the detection angle is analogous to 
performing a series of conventional DWS measurements 
using a series of sample cells with varying thickness or 


with varying tracer particle concentrations. 

We have further illustrated the usefulness of our 
method for DWS-based microrheology on a soft solid, 
gelatin, for which we obtain a very good agreement with 
macroscopic oscillatory rheology experiments. Moreover, 
data recorded for different detection angles enable an im¬ 
portant consistency check for the microrheology measure¬ 
ments and our results illustrate that the accuracy of such 
microrheology measurements can be improved by aver¬ 
aging over measurements obtained at different detection 
angles. 
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Error analysis for DWS and microrheology data 

In the following we provide a brief analysis aimed at estimating the relative experimental errors associated with 
our DWS-based microrheology measurements. 

In DWS, the dynamics of tracer particles is quantified in terms of temporal autocorrelation functions, from which 
the time-dependent mean-square displacement hereafter abbreviated as MSD(f) ) of the particles can 

be calculated. Using a generalized Stokes-Einstein relation, the MSD’s are then converted to viscoelastic moduli, 
with the magnitude of the complex modulus |G*(a;)| approximately proportional to the inverse of the mean-square 
displacement at a time scale t = 1/uj. 

The relative error in the modulus, AG*/G*, is therefore in good approximation the same as the relative error in 
the mean-square displacement AMSD/MSD. We identify two main sources of error that ultimately determine this 
relative error in the measured mechanical response: 

1. Statistical errors in the intensity correlation function as a result of the finite measurement duration. 

2. For nonergodic samples, additional errors are introduced via the Pusey averaging procedure used to estimate 
ensemble-averaged temporal autocorrelation functions. These errors are introduced via the relative error in the 
intensity ratio Y = between the time-averaged and the ensemble-averaged scattering intensities. 


a. Statistical errors in the intensity correlation function 


Given a measurement duration T, the statistical error in the intensity correlation function g 2 {t) at lag time t can 
be expressed as m 


Ag2{t) « 6 



(9) 


Given the Siegert relation gi{t) = yj 52 (t) — 1, this translates to an error in the field correlation function gi{t) as 
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51W 


( 10 ) 


For the purposes of this error analysis, we assume a simplified relationship between the field correlation function 
gi(t) and the mean-square displacement MSD(t), gi{t) « e“*^o/3itMSD(t)^ with ko the wave vector of the laser light, s 
the average path length of photons, and I* the transport mean free path. 

31* 1 ^ 

91 


We can then express the mean square displacement as MSD(<) « 31* \n{gi{t)/ (fco^s)- With 9MSD/9gi 
li 


91 


MSD/( 5 i ln( 5 i)) and A^i ^ we can now express the total relative error in the complex modulus as 


AG*/G* « AMSD/MSD : 


1 (9MSD 


A^i 


1 


( 11 ) 


MSB dgi gilii{gi) gi\ T 

b. Errors introduced during the Pusey averaging procedure 

For nonergodic samples, additional errors are introduced via the Pusey averaging procedure used to estimate 
ensemble-averaged temporal autocorrelation functions. These errors are introduced via the relative error in the 
intensity ratio Y = ^ between the time-averaged and the ensemble-averaged scattering intensities. The resulting 
error in the mean square displacement is 


AMSD^^^^AY 

ogi dY 


( 12 ) 


where , as derived above. 

agi 9 iln( 9 i)( 

To arrive at an expression for the second term in Eq |12[ we write down the relationship between the field correlation 
function gi{t) and the intensity ratio Y as explained in the main manuscript. 


, ^ Y -1 1 r . ^ 91 

51 (1) = + y [52(0 - cr ] 


(13) 


where 32 (i) = 1 + —- is the time-averaged intensity autocorrelation function normalized by the coherence factor 

/3, and = ^ 2(^)1 characterizes the short-time intercept of 52 (i)- This yields the second term in Eq 
1 - \/ff2 - 


1 

Y2 


12 


as 


dgi 

dY 


and with \/g 2 — ~ ^ YYgi — Y, this results in 

1-51 


^51 

dY 


Y 


(14) 


The total relative error introduced by the intensity ratio Y is thus 

1 -5i 


AMSD 


5iln(5i 


AY 


(15) 


Plotting the function f{x) = we observe |/(a:)| « 1 for values of x close to 1. For our solid-like samples. 


gi is close to 1 at long times, and thus 


AMSD ^ AM 
MSD ~ Y ■ 


I- 9 I 


9 lln( 9 l) 


1; therefore in this case we can make the simple approximation 


c. Total estimated error 


The total estimated error as a result of the above two main sources of error can then be written as 

3 



1-51 

AY 

V r ^ 

51 ln(5i) 

Y 


AG*/G* « AMSD/MSD 


(16) 
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Of these two contributions to the experimental error, the latter usually dominates, provided that the measurement 
duration is sufficiently long, such that the factor becomes small enough. Nevertheless, for gi{t) sufficiently close 
to 1, the term would become very large, and eventually lead to the first term becoming dominant. 

Indeed, we expect the errors to be highest for cases where gi (t) is very close to zero (as a result of the factor ^^2 in(gi) 
or , respectively), or very close to unity (as a result of the factor 


1. Code for data analysis of DWS in a cylindrical cell 


a. Full codes (written in Matlab) for analyzing DWS data in a cylindrical cell 

The complete numerical codes used for analyzing dynamic light scattering data using the approach outlined in the 
manuscript can be obtained on the author’s website <www.mate.tue.nl/^wyss> or on request by sending an email 
to Hans Wyss at H.M.Wyss@tue.nl, Please also address any questions regarding the code and/or data analysis to the 
same email address. 

The specific code for the random walk simulation for calculating the path length distribution P{s) is listed below. 


b. Code for calculating the path length distributions 

The main code for calculating the path length distribution P(s) is a simple random walk simulation with step length 
I*, as described in the main manuscript. Below is the C-code (and the MEX function called by our main Matlab 
code) that accomplishes this task; the results are kept track of in the output array y, which for each bin corresponding 
to a segment of detection angle and a segment of path length keeps track of the number of photons exiting the cell 
within the corresponding angle range and within the corresponding path length range. Angular bins evenly divide 
the angular space between 0 and 180 degrees; we usually choose 180 bins of 1 degree width. Path length bins are 
also linearly spaced, with 300 bins total and the 100*’*' bin corresponding to a path length of (L(0)/P)^ in units of Z*, 
where L(9) is the distance between the entry point and the exit point of the simulated photon, as a function of the 
detection angle 0. 


Listing of “pathlengthsCyl_P_s.c”: 

1 ^include ’’mex.h” 

2 ^include <stdio.h> 

3 ^include <stdlib.h> 

4 ^include <time.h> 

5 ^include <math . h> 

6 

7 A 

8 * 

9 * Function that calculates the pathlength distribution for diffusion of photons 

10 * through a eyUnder 

11 * 

12 */ 

13 

14 void pathlengthsCyl ( double y[] , double x[] , size_t mrows, size_t ncols) 

15 { _ ^ 

16 int ii ,m, n,step,anglebin, Nbins , N_iter , index 1 , index2 , index 3 ; 

17 double xx , yy , zz , dxx , dyy , dzz , 1st ar ,LL, pi , deltaz , weight , dist ; 

18 

19 lstar=x [0] ; 

20 Nbins=x [ 1 ] ; 

21 N _it er=x [ 2 ] ; 

22 pi=3.14159265359; 

23 








14 


24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 


40 

41 


42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 


//srand( (unsigned) time ( NULL ) ); 
srand ( rand() " time (NULL) ); 

for (m=0;m<mrows ;in-|-+) 

{ 

for (n=0;n<ncols ; n++) 

{ 

indexl=(m% mrows)+mrows*n; 
y [ indexl ] = 0; 

} 

} 

for (n = 0;n<ncols ; n++) 

{ 

dist=sqrt ( 

(—1— cos(pi/Nbins*n))*(—1—cos(pi/Nbins*n))+sin(pi /Nbins*n) *sin(pi/Nbins*n) 

); 

indexl=2+mrows*n; 

y [ indexl]= cei 1 ( dist * dist / 1st ar / Ist ar /100) ; //write column 2: the width of 
each bin (set as one hundreds of the expected average number of 
scattering events.) 

} 


for( ii=0;ii<N_iter ; ii ++) 

{ 

xx=—1+1 s t a r ; 
yy = 0;zz=0; 
step = 1; 

while (( xx*xx+yy*yy) <1) 

{ 


// Random unit vector of length Istar: 
dxx = 2; dyy = 2; dzz = 2; 

while (dxx*dxx+dyy*dyy+dzz*dzz >1) // make sure (dxx , dyy , dzz) is a 
vector of random direction . 

{ 

dxx=2*(double) (rand 0 ) / (double) (RANDJMAX)—1; 
dyy=2*(double) (rand 0 ) / (double) (RANDJMAX)—1; 
dzz=2*(double) (rand 0 ) / (double) (RANDJMAX)—1; 

} 

LL=sqrt(dxx*dxx+dyy*dyy+dzz*dzz); 

dxx=dxx/LL*Istar;dyy=dyy/LL*Istar;dzz=dzz/LL*Istar; 

// Propagate by a step Istar: 
xx=xx+dxx ; yy=yy+dyy ; zz=zz+dzz ; 
step++; 

} 

LL=sqrt(xx*xx+yy*yy); 
xx=xx/LL;yy=yy/LL;zz=zz/LL; 
anglebin=cei 1 ( acos (xx) *(Nbins —1)/pi) ; 

indexl=mrows* anglebin ; 
deltaz=pi /2/Nbins ; 

weight=er f (3* delt az / Is t ar / sqrt ( step ) ) ; //account for probability in 
z—direction to hit area around zz=0 
if (weight==0) mexPrintf (” erf : %f\n” , weight) ; 
y[indexl] = y[indexl] + weight ; 
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76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 } 

103 

104 

105 

106 

107 

108 void mexFunction ( int nlhs , mxArray *plhs [] , 


109 

int nrhs , const mxArray *prhs [] ) 


no { 



111 

double *x,*y; 


112 

double 1st ar ; 


113 

int Nbins ; 


114 

size_t mrows , ncols , mrows_ont , ncols_out ; 


115 



116 

/* Check for proper number of arguments. */ 


117 

if ( nrhs ! = 1) { 


118 

mexErrMsgIdAndTxt ( ’’MATLAB: timestwo : InvalidNum 

Inputs” , 

119 

’’One input required.”); 


120 

} else if ( nlhs >1) { 


121 

mexErrMsgIdAndTxt ( ’’MATLAB: timestwo : maxlhs” , 


122 

’’Too many output arguments.”); 


123 

} 


124 



125 

/* The input must be noneomplex doubles.*/ 


126 

mrows = mxGetM(prhs [0]) ; 


127 

ncols = mxGetN) prhs [0]) ; 


128 

if( ! mxIsDouble ( prhs [ 01) I I mxIsComplex ( prhs [0]) ) 

{ 

129 

mexErrMsgIdAndTxt ( ’’MATLAB: timestwo : inputNotRealScalarDouble” 

130 

’’Input must be noncomplex doubles.”); 


131 

} 


132 




// Update average number of steps and corresponding added weights: 
index2 = (l % mrows)+mrows* ang 1 ebin ; 
y[index2] = y[index2] + step*weight; 

index3=floor ( step /y [ anglebin *mrows + 2] )+3+mrows* anglebin ; 
if ( floor (step/y[anglebin* mrows + 2 ]) +3<mrows) 

{ 

y [ index3] = y [ index3] +weight ; // add value to bin 

} 

} 

for ( anglebin =0; anglebin <Nbins ; anglebin++) 

{ 

indexl=mrows* anglebin ; 

index2 = (l % mrows)+mrows* anglebin ; 

if(y[indexl]==0) 

{ 

// mexPrintf (” y [ indexl ] is zero. indexl:%i yfindexl]: %f, y[index2]: 
%f\n”,indexl, y[indexl], y[index 2]); 

} 

if (y [ index2] >0) 

{ 

y[index2]=y[index2]/y[indexl ] ; 

} 

} 
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133 /* Assign pointer for the input matrix.*/ 

134 X = mxGetPr ( prhs [ 0 ]) ; 

135 lstar=x [0] ; 

136 //x f 1 ] = flo or (Nbins) ; 

137 Nbins=x[l]; 

138 

139 /* Create matrix for the return argument. */ 

140 

141 mrows_out =303; 

142 // row 0: average intensity for this angle bin; 

143 // row 1: average path length for this angle bin. 

144 // row 2: bin width for this angle bin. 

145 // rows 3—302: path length distribution for this angle bin. 

146 ncols_out=Nbins; /* one column for each angle segment. The columns give a 

counter for each */ 

147 plhs [0] = mxCreateDoubleMatrix (( mwSize) mrows.out , (mwSize) ncols_out , mxREAL) ; 

148 

149 /* Assign pointer for the output matrix. */ 

150 

151 y = mxGetPr( plhs [0]) ; 

152 

153 /* Call the pathlengths subroutine . */ 

154 pathlengthsCyl(y,x, mrows_out ,ncols_out); 

155 } 

pathlengthsCyLP_s.c 
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